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Experiments that discuss influence of noise to H. Seidel and H. Herzel dynamics model of human 
cardiovascular system are presented. Noise is introduced by considering stochastic delays in response 
to the sympathetic system. It appears that in the presence of the noise 10 s heart rate oscillations 
connected with Mayer waves are preserved. Moreover the heart rate becomes approximately nor- 
mally distributed (even in unstable phase of original Seidel Herzel model) , similarly like the real RR 
intervals data. 
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I. INTRODUCTION 

The human cardiovascular system is controlled by sev- 
eral neural and hormonal mechanisms. Among these 
mechanisms the autonomic nervous system via barore- 
ceptors activity provides blood pressure control 1] . This 
control is realized by two parts of the autonomic sys- 
tem acting antagonistically: the sympathetic subsys- 
tem which activity accelerates the heart rate and the 
parasympathetic subsystem which activity elongates the 
interbeat intervals. The finite propagation speed of neu- 
ral impulses induce several time delays. Time delays are 
a very common feature of physiological systems Q • Their 
incorporation is essential for a realistic modelling. 

Modelling of the baroreceptors-cardiac reflex dynamics 
has a long history, see e.g. 0, 0, 013 ■ 0ne of tne 
still widely investigated models was proposed by DeBoer 
et al. 3]. In this model blood pressures, respiration, 
peripheral resistance and cardiac interbeat intervals are 
represented by a simple beat-to-beat relations. Ottesen 
4] modelled chronotropic (heart rate regulation) and in- 
otropic (contractility of ventricle regulation) parts of the 
baroreflex-feedback mechanism. This proposition based 
on expanded Windkessel model, describes both sympa- 
thetic and parasympathetic nervous system, with one 
sympathetic time delay. The system can switch between 
being stable and oscillatory, and vice versa, with changes 
of sympathetic time delay. Eyal et al. [5( studied sim- 
plified linearized version of DeBoer model. They found 
that increase of sympathetic gain or decrease of vagal 
gain parameters leads to Hopf bifurcation and sustained 
oscillations. 

The idea of Seidel and Herzel proposition, called SH 
model, consists in entering the internal dynamics of one 
heartbeat 0. It is achieved by introducing a phase of 
a sinus node (the first pacemaker of the heart) and then 
considering the sensitivity of the sinus node phase change 
to autonomic neurotransmitter kinetics and vascular dy- 
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namics. Moreover, the SH model is strongly physiologi- 
cally grounded. Parameters of the model take experimen- 
tally known values and the solution restores the rhythms 
of heartbeat with a period of about 1 second, respiration 
with a period of 4 seconds and Mayer waves with a period 
about 10 seconds. Therefore this model is often used in 
modelling other phenomena related to heart dynamics, 
ejj., synchronization between respiration and heart rate 
@ or presence of 1/f fluctuations in interbeat interval 
signals 0. 

Since interactions in SH model involve nonlinear de- 
pendencies then chaotic aspects of dynamics are present. 
An increase of the delays in conveying of sympathetic ner- 
vous system signals leads via Hopf bifurcation to heart 
rate oscillations [13J . Such altered dynamics as the result 
of bifurcation is present in many mathematical models 
of physiological systems. Especially systems with feed- 
backs and time delays can exhibit both oscillatory and 
chaotic dynamics |l2j. Bifurcation with time delays as 
bifurcation parameters is observed also in other models 
of baroreflex 0HG3. 

It is known that cardiovascular signals behave deter- 
ministically, but noises of various origin are also present 
in the system [Tpj. For example, the basic neurotransmit- 
ter of the sympathetic nervous system - noradrenaline, is 
distributed to the heart and vessels via diffusion. There- 
fore this process should be represented by slow stochas- 
tically dependent dynamics. It would be interesting to 
change SH model by considering stochastically driven pa- 
rameters describing delays in signal transfers. The goal 
of our investigation is to examine the stability of the solu- 
tion of SH model when the delay parameters are stochas- 
tic variables. 

We focus our interest on the bifurcation regime. This 
regime demands careful dealing with model parameters. 
We have found that there is a parameter for which the 
solution of the original SH model has no physiological 
meaning. This parameter is involved in the construction 
of the sigmoidal function which is responsible for Wind- 
kessel time constant dependence on large vascular con- 
centration of noradrenaline. To keep the solution within 
physiologically reasonable limits we propose to change a 
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value of this parameter. The change proposed by us does 
not influence the model solution in the regular regime. 
However it allows to find solutions which are numerically 
stable and physiologically rational for wide range of time 
delays. 

For the reasons described above we present the ba- 
sic model and our modification in details, Section 2, 
and then discuss thoroughly obtained solution, Section 3. 
Section 4 contains simulations when delays are stochas- 
tically changed. 



II. SH MODEL 



Sympathetic influence on the phase velocity of the 
sinus node: 



fs = l + kf a [ (c cNa + 

(c c Na —CcNa) 



-cNa 



nsd(6) 



(k^ Na = 1.6,c cNa = 2.0,n cNa = 2.0 - parameters) 

Parasympathetic influence on the phase velocity of 
the sinus node: 



SH model provides relations between different mecha- 
nisms involved in the baroreceptor-cardiac reflex in the 
set of the following equations: (see FigQ] for the visual- 
ization of relations between variables) 

- Baroreceptor activity v : 



is b = k 1 (p-pW) + k 



dp 

a 



(i) 



(p - blood pressure; k\ = 0.02mmHg 1 ,k2 = 
0.00125s mmHg -1 ,p(°) = 50mmHg - parameters) 

— Sympathetic activity v s : 

v s = max(0, i/W - k\v b + k r s \ sm(%f r t + A#)|) (2) 

(fr — 0.2s -1 - respiratory frequency; isi ^ = 
0.8, k b s = 0.7, k r s = 0.1, A(f> r s = 0.0 - parameters) 

— Parasympathetic activity v v : 

u p = max(0, + k\v h + k r p \ wa.(nf r t + A^)|) (3) 

(f r - respiratory frequency; — 0.0, kp — 

0.3, kp — 0.1, Acfrp — 0.0 - parameters) 

- Cardiac concentration of sympathetic transmitter 

CcNa'- 



dc C N a CcNa 



dt 



TcNa 



+ kl Na V s {t-0cNa) (4) 



(T C Na — 2.0s, kc cNa — 1.2 - parameters, C Na - time 
delay) 



Sinus node phase ip: 



dip 1 



(5) 



(fs,fp - sympathetic and parasympathetic influ- 
ences; r(°) = 1.1s - parameter) 



fp = !- k l{ v p , { 



(Vp -VpM p ) ^.n p ' "n,, \F((P) (7) 



( v pfi v = v p(t — 6p) - delayed parasympathetic activ- 
ity, F(ip) - phase effectiveness curve, kR = 5.8, v p = 
2.5, n p = 2.0- parameters) 

Phase effectiveness curve F: 

Cardiac contractility without saturation S'f. 

S[ = + k C s C cN a + 4^-1 (9) 

(Tj_i - duration of previous heart period; — 
25mmHg, kg = AQmmHg, k s = lOmmHgs -1 - pa- 
rameters) 



Cardiac contractility with saturation Sf. 

o'ns 

S i = S' i + (S-Sl). 



(S = 70mmHg,ns = 2.5 - parameters) 
Windkessel time 'constant' t v : 

T v = - T~[ C vNa + 



(10) 



J vNa 



(c v Na C v Na)-^n^ N ,, , v 

C vNa ">" C vNa 



(11) 



(Tv = 2.2s,Tv = 1.2s, C v Na = 10.0, n v Na = 1.5 - 
parameters) 

Model modification: Fig. [21presents the dependence 
of t u - Windkessel time 'constant' on c V Na, follow- 
ing Eq. ^] Such dependence admits t v to be neg- 
ative. However since t v is responsible for the blood 
pressure decay in aorta during diastolic part of the 
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FIG. 1: The scheme of SH model of the baroreceptor control loop. Arrows depict relations between variables. 




FIG. 2: The dependence t v on c V Na- tI = t v (c v no) for 

CvNa = 10.0,7^ = r v (c vNa ) for C vNa = 1.0 



heartbeat cycle, see eo.lfTljl. hence it must be posi- 
tive. Therefore we propose the following change in 
the sigmoidal function parameter: let c V Na = 1.0 
in place of c V Na = 10.0. This change becomes im- 
portant only in a case when delays 9 C N a , 9 C Na are 
greater than 3.0. 



Vascular concentration of sympathetic transmitter 

CvNa'- 



dc v N a C v Na 



(It 



TvNa 



+ K vN M t ~ vNa) (12) 



{TvNa = 2.0s, k s CvNa = 1-2 - parameters, V Na - time 
delay) 



Blood pressure during the systolic part of the heart 
cycle: 



p = di-i + Si -exp < 1 



' sys 



* sys 



(13) 



(<ii_i - diastolic pressure during the systolic part of 
the heart cycle, ti - time of last contraction onset, 
i~sys = 0.125s - parameter) 

Blood pressure during the diastolic part of the 
heart cycle: 



dp 
~dt 



P 



At) 



(14) 



4 



III. NUMERICAL SIMULATION OF THE 
MODEL 

The adapted Runge-Kutta method of fourth order with 
a constant step size (h=0.001s) was used. The ring 
buffers are introduced to store history of sympathetic and 
vagal activities. The first 500 seconds of evolution are ne- 
glected, as the transient time, and then the results are 
collected and analyzed. The results with parameters set- 
ting as in the original Seidel-Herzel proposition are pre- 
sented in Fig. |31 Subsequently, starting from the top of 
Fig. |3| we show: 

• The blood pressure p which is characterized by 
rapid increase during the systole, and gentle de- 
crease during the diastole. Minimum (diastolic - 
about 80 mmHg) and maximum (systolic - about 
140 mmHg) values of this wave nicely reproduce 
physiological one's. 

• Baroreceptors activity Vb that basically is propor- 
tional to the blood pressure. Hence the waves in 
the first and the second panels are similar to each 
other. 

• The parasympathetic activity v p that affects the 
heart with the time delay 9 P and because of this 
delay the oscillations of the activity are shifted in 
time. Since transmission of the signal goes fast 
(mechanism of transmission is different from dif- 
fusion) then the delay can be assumed as constant 
and relatively small, i.e., 9 P = 0.5s. 

• The sympathetic system activity v s enters into the 
model with two delays: 9 C Na and 9 V Na correspond- 
ing to its influence to the heart and vascular system, 
respectively. Delays 9 C Na, @vNa represent times 
needed by transmitters to transfer the informa- 
tion by diffusion. The activity affecting the heart 
v s (t — 9 c na) is presented in the figure. 

• Cardiac c C Na and vascular c V Na noradrenaline con- 
centrations which take the same values. 

• The phase of the sinus node tp which mimics heart 
beating (integrate-and-fire-modcl) and establishes 
the rhythm of the whole system. All variables of 
the model oscillate with a period about 0.8 s. 

• The sequence of heart periods Ti that represents the 
time interval between two consecutive heart beats. 

Our results are slightly different from results presented 
in . To explain let us observe that there is some incon- 
sistency in the results of 0. For example, it is easy to 
estimate (directly from equations J2J) and ©) that the 
maximum and minimum values of v p and v s plotted in 
Fig. 4 of should be different: v p should be smaller 
and maxima of v s should be greater than it is shown. 

Sympathetic and parasympathetic activities are influ- 
enced by respiratory neurons. This influence is described 



by the sine function in Eq. (0) , J3J . Blood pressure and 
heart period variations caused by respiratory are called 
respiratory sinus arrhythmia (RSA). Respiratory modu- 
lation of neural activities is clearly visible in Fig[3Ia). All 
waveforms are modulated with the frequency of respira- 
tion. As in our further simulations we concentrate on 
mechanisms other than respiration then the sine func- 
tion is replaced by the average value, and the respiratory 
modulation of neural activities is neglected. Fig. |3Ib) 
shows the solution of the model with the neglected res- 
piratory modulation. In the absence of RSA the heart 
period is approximately constant but it is not settled to 
a one point. Let us call this solution a regular heart rate. 

We present two representative examples to illustrate 
how time delays influence the solution of the model. The 
first bifurcation diagram is obtained for 9 v ^ a = 1.65s 
and increasing 9 C M a , see Fig. and the second one is 
for 9 V N a = 3s and increasing # c at q , Fig|S] One can easily 
separate regions of the regular heart rate (9 V Na = 1-65: 
9 C N a < 2s and 9 V Na = 3: 9 C Na < 0.6s) from the heart 
rate with oscillations {9 v ^ a = 1.65: 2 < 9 c ^ a < 3.5s and 
9 V Na = 3: 0.6 < 9 C N a < 3.5s). The influence of delays 
is not symmetric: the setting 9 V Na = 1-65, 9 C N a = 3.0 
leads to heart periods 0.45 < T < 1.1 while the setting 
9 vNa = 3., 9 cNa = 1.65 leads to 0.55 < T < 0.9. 

IV. STOCHASTIC APPROACH 

Let us assume that a time, needed for a signal to trans- 
fer information from the brain through the sympathetic 
nerves to the heart, varies stochastically around some 
mean value. To incorporate this assumption let us in- 
troduce random variables in place of fixed delays in the 
activity function of the sympathetic nervous system. So 
Eq. (0J and l|12|l obtain the following form: 

(1 f'M C^f. 

= — + klu s {t-{9 sr +^)) (15) 
at t* 

where are random variables with uniform distribution 
and G [— £*,£*], and * represents either cNa cardiac 
or vNa vascular noradrenaline concentration. The ran- 
dom numbers are generated by standard 48-bit arith- 
metic functions accessible in C compilers in Unix systems. 
Noradrenaline needs to be removed from the synaptic gap 
and broken down at another location. This process takes 
much longer time than a time step of model simulation. 
Noradrenaline concentration acts like a buffer so the dis- 
tribution of sympathetic impulses over a heart cycle is 
not very important |13| . Therefore we assign random val- 
ues to delays for each cardiac cycle separately. One can 
model the noise by a continuous variable also. In this 
case, at each numerical step a different random value to 
the delay is set. The randomness put into the system in 
such a way does not cause significant change in solutions 
of the model. 

To avoid effects related to respiration we again consider 
the system with constant respiration only. The influence 
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FIG. 3: Time series of the SH model variables. Time delays are fixed O c jva=O t ,]Va = 1.65s. (a). The influence of respiration can 
be seen in modulations of waveforms (RSA) . (b) . Time series when the respiratory modulation of naural activities is neglected 
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FIG. 4: The delay in transmitting signals from the brain to 
the heart 8 C Na is increasing, while the delay between the brain 
and the blood vessels is constant 8 V Na = 1.65s 



FIG. 5: The delay in transmitting signals from the brain to 
the heart 9 C Na is increasing, while the delay between the brain 
and the blood vessels is constant 8 V Na = 3s 



of the stochastic perturbation to the system will be ex- 
amined when the heart beats regularly: 6 C Na — &vNa = 
1.65s, and when the heart rate oscillates: 9 C Na — &vNa = 



3s. Statistical properties of the heart rate will be pre- 
sented by histograms - the probability that a given period 
is present in a series. For comparison, let us show few 
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FIG. 6: Histograms of time series of RR intervals from healthy 
subjects [l5| . To study the similarity to normal distribution 
the log scale is applied. 
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examples of histograms made from real heart rate series 
(so-called normal sinus rhythm series), Fig. |SJ Although 
each of the histogram is different from others but the 
parabola-like shape in log-scale suggests that the Gaus- 
sian distribution of heart rate intervals can be used as 
the first approximation to the histograms. 



A. Stochastic perturbation to regular heart rating 

If the stochastic time delays are considered then the 
sympathetic activity becomes irregular what affects oth- 
ers variables — the heart rate also. Fig. [7| shows the 
time series obtained when the stochastic delays are ap- 
plied. The regular waveforms received in deterministic 
case for the sympathetic activity, compare to Fig. |3|d, 
now, at the presence of stochastic delay, are varying ir- 
regularly. The unperturbed evolution provides the his- 
togram with one value of heart period only, see Fig. [S] 
Depending on the magnitude of added noise this line is 
broadened. Notice, that at high level of the noise, namely 
if 0.65s < C N a ,9 v pf a < 2.65s, the histogram shape re- 
minds distributions of real series. 



B. Stochastic perturbation to irregular heart rating 

Fig. & presents unperturbed time series for 9 c ^ a = 
d V Na = 3 while Fig. |5Jd shows series when delays are 
2 < C N a , OvNa < 4 taken at random. 

The SH model is able to generate sustained 10 s 
rhythms in the bifurcation regime [t| which can be re- 
lated to Mayer waves known from physiology. However 
the statistical properties of heart period are strange. The 
histogram of the heart beat intervals in case of the unper- 
turbed dynamics takes the form of a double peaked distri- 
bution with sharp decay of the left and right wings, Fig. 
ITU1 The shortest and the longest periods appear with 
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FIG. 7: Time series obtained by simulating the model with 
stochastic delays - regular heart rate case. 
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FIG. 8: Histograms of heart period T - regular heart rate case. 
The unperturbed evolution (top-left panel) and evolution with 
increasing level of noise. To study the similarity to normal 
distribution the log scale is applied. 



highest probability. The heart period looks like switch- 
ing between these two values. Such a distribution is not 
resembling any distribution found for real RR time series. 

Stochastic noise influences the heart rate but the main 
oscillations of the deterministic solution such as Mayer 
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FIG. 9: Time series obtained by simulating the model (a) without and (b) with stochastic delays - oscillating heart rate case 
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FIG. 10: Histograms of heart period T - irregular heart rate 
case. The unperturbed evolution (top-left panel) and evolu- 
tion with increasing level of noise. To study the similarity to 
normal distribution the log scale is applied. 



waves are partially preserved, see Fig|^b). However, the 
distribution of the intervals is wider — the intervals from 
peaks are less probable. When the noise level is such that 



conditions for a regular heart rating emerge, i.e., 1.5 < 
6cNa,8vNa < 4.5, then the distribution almost restores 
the normal shape, see Fig. ITU1 right-bottom plot. 



V. CONCLUSION 

In this paper we investigated the dynamic properties of 
stochastic nonlinear model of the cardiovascular system. 
The stochastic model was obtained by allowing random 
changes to bifurcation parameters of the Seidel-Herzel 
model. Using the computer simulations the role of the 
stochastic perturbations was studied. It is shown, that 
stochastic delays lead to time series with probabilistic 
properties similar to those obtained from real time series 
for healthy humans. Our results ensure that incorporat- 
ing with random variables and in consequence consider- 
ing stochastic differential equations in place of determin- 
istic ones, we gain more realistic description of the car- 
diovascular system. Moreover, the introduced change in 
sigmoidal function parameter allows to explore the model 
for wide range of time delays. As SH model is widely used 
to reproduce key features of human heartbeat dynamics 
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(e.g. |E 0) the possibility of carrying out simulation for /5400-5-0166-5 

wide range of parameters is very important for further 

study. 
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